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Time-distance helioseismology and related techniques show great promise for probing the structure and dynamics of the 
subphotospheric layers of the Sun. Indeed time-distance helioseismology has already been applied to make inferences 
about structures and flows under sunspots and active regions, to map long-lived convective flow patterns, and so on. Yet 
certainly there are still many inadequacies in the current approaches and, as the data get better and the questions we seek 
to address get more subtle, methods that were previously regarded as adequate are no longer acceptable. Here we give a 
short and partial description of outstanding problems in local helioseismology, using time-distance helioseismology as a 
guiding example. 
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T 1 ' 1 Introduction 
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Local helioseismology has grown spectacularly in the past 
decade. Global helioseismology, using principally the fre- 
quencies of global modes, has enabled us to image the radial 
variation of the hydrostatic structure through most of the 
solar interior, and similarly the radial and latitudinal vari- 
ation of the rotation rate. But some interesting features of 
the solar activity cycle - sunspots, active regions - are by 
definition horizontally localized and influence global-mode 
frequencies only very subtly. Global modes are not the best 
instruments to probe these features. Local helioseismology 
is much better suited to this task. A number of local helio- 
seismological techniques and approaches have been devel- 
oped. These include time-distance helioseismology (Duvall 
et al. 1993), ring analysis (Hill 1988; Antia & Basu 2007), 
acoustic holography (Lindsey & Braun 2000) and statistical 
waveform analysis (Woodard 2002). In this paper we are 
concerned mostly with time-distance helioseismology: for 
a recent review, see Gizon & Birch (2005). 

Time-distance helioseismology proceeds by cross- 
correlating measurements (typically Doppler velocities) of 
the oscillations at different locations on the solar surface, 
fitting the cross-correlation to determine a 'travel time' 
for waves propagating between those different locations, 
and then interpreting those travel times via some inversion 
procedure to infer something about the subsurface condi- 
tions affecting the wave propagation. What issues might 
be addressed using helioseismology regarding, for exam- 
ple, sunspots? The list includes the structure of a sunspot 
both horizontally and with depth, the flows around the 
sunspot and how these vary with depth and distance from 
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the sunspot. How can we decide what is the best technique 
to employ to exploit the data at hand? Time-distance helio- 
seismology has proved attractive for a number of reasons. 
It is quite intuitive. It is fairly robust, in that the main mea- 
surement is a fitted travel time rather than more detailed as- 
pects of the cross-correlation function or wave field which 
could be more sensitive to the manner of mode excitation. 
It achieves a much better resolution than has hitherto been 
achieved by ring analysis, a technique which requires ob- 
servations taken over a relatively large patch of Sun in order 
to get adequate wavenumber resolution. In contrast, time- 
distance spatial resolution may get down to the wavelength 
of the waves being used in the analysis. 

Woodard's statistical waveform analysis has been hailed 
as using more completely than any of the other local he- 
lioseismic methods the information in the observed wave- 
field, and therefore having the potential to provide more de- 
tail about the underlying subphotospheric conditions. Yet at 
present the technique is certainly not well developed, there 
is the issue of whether it would require impractical amounts 
of computational resource and how much additional infor- 
mation can be extracted compared to other methods. 

2 Dopplergrams 

The main observational signal used in helioseismology is 
the Doppler velocity, which is obtained by combining inten- 
sity measurements in the core and wings of a photospheric 
absorption line (e.g. Toussaint, Harvey & Hubbard 1995; 
Scherrer et al. 1995). Solar intensity images, or filtergrams, 
are recorded with a CCD camera over short integration in- 
tervals to provide one Doppler velocity image every minute. 
The study of solar oscillations requires long time series of 
Dopplergrams - sampled images of the solar line-of-sight 
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velocity convolved by the point-spread function (PSF) of 
the telescope. 

Several instrumental effects complicate the inter- 
pretation of Doppler velocities and lead to system- 
atic errors in the seismic measurements. Understand- 
ing these effects is an important and challenging task. 
|Korzennik, Rabello-Soares & Schou~(2 004) have discussed 
imperfections in the MDI image plate scale, the im- 
age distorsion, and temporal changes in the optical 
system. The MDI PSF, which has been estimated by 
Tarb ell, Acton & Frank (1996)| near disk center, is not ax- 
isymmetric, implying a different sensitivity to waves prop- 
agating in different directions. Woodard et al. (2001) stress 
that a good knowledge of the PSF is critical in order to be 
able to infer by deconvolution the true energy of solar oscil- 
lations versus horizontal wavenumber. 

Systematic variations across the field of view are, of 
course, not only due to instrumental effects but also to 
purely geometrical effects including foreshortening - the re- 
duced sensitivity to short wavelengths away from disk cen- 
ter in the center-to-limb direction - and the projection of 
wave motion (mostly vertical) onto the line of sight. In ad- 
dition, an error in the orientation of the image will affect 
the mapping into the heliographic coordinate system. In- 
deed Giles (1999) showed that a 0.1° P-angle error can lead 
to a measurable ~ 5 m/s spurious flow in the north-south 
direction caused by a leakage of the 2 km/s solar rotation. 
|Giles (1999)| also found a one-year periodicity in the MDI 
travel-time measurements, which he attributed to an error 
in the determination of the direction of the solar rotation 
axis. The high-precision alignment and merging of Doppler- 
grams from the GONG ground-based network of telescopes 
is discussed by |Toner (200 1)[ 

A good calibration of the velocity signal requires 
a physical model of the absorption line profile from 
which it is derived. Structural inhomogeneities in the 
solar atmosphere, especially magnetic field concentra- 
tions, affect the shape and formation height of a spectral 



line. |Wachter, Schou & Sank arasubramanian (2006) and 
Raj aguru et al. (2006b) have proposed corrections to the 



standard calibration to account for systematic errors due to 
spectral line shape changes in solar active regions. These 
studies are important in order to make a distinction between 
a biased estimate of the Doppler signal and a physical per- 
turbation to the seismic wave field. A source of concern is 
the possibility of a cross-talk between the measured veloc- 
ity and magnetic signals. Besides the magnetic field, tur- 
bulent flows also affect line profiles. A radiative transfer 
calculation of spectral line profiles formed at two different 
heights in the solar atmosphere is presented in this volume 
by Haberr eiter, Finsterle & Jefferies (2 007) 

The observational systematic effects mentioned above 
have a direct influence on the interpretation of local helio- 
seismological measurements. Correcting for all these effects 
is a formidable, but highly desirable enterprise. 



3 Travel times 

In local helioseismology, the governing equations of so- 
lar oscillations are often written in plane-parallel geometry. 
This approximation is typically valid over a small fraction 
of the solar disk for horizontal wavelengths that are much 
smaller than the solar radius. Ignoring most of the system- 
atic errors discussed in the previous section, we may write 
the Doppler signal used in helioseismology as the line-of- 
sight component of the wave velocity, 



(x,t) =PSF*[e-dtt(x,z ,t)} 



(1) 



where a; is a horizontal position vector on the solar surface 
(defined by height Zq), t is time, i is a unit vector pointing in 
the direction of the observer, and £ is the wave displacement 
vector. Convolution by the PSF is included. It is customary 
to consider the data cube <fi(x, t) in a reference frame which 
is co-rotating with the Sun to remove the main component 
of rotation. 

In time-distance helioseismology, the oscillation signal 
is filtered in Fourier space to obtain 



i>(k,L>)=F(k,L>)4>(k,U>), 



(2) 



where k is the horizontal wave vector, ui is the angular fre- 
quency, and F is a (real) filter function chosen by the ob- 
server to remove granulation noise and to select parts of the 
wave propagation diagram. A standard procedure is to use 
phase-speed filters |(D~u vail et al. 19 97) This choice is moti- 
vated by the fact that acoustic waves with the same horizon- 
tal phase speed w/k are confined to the same acoustic cav- 
ity |(Bogdan 1997)| The notion of optimal filtering remains 
to be defined, but is appealing |(Couvidat & Birch' 2006) 

The fundamental time-distance computation 
|(Duvall et al. 1993)1 is me cross-covariance function 
between the filtered Doppler velocities at two surface 
locations X\ and a; 2, 



C(x 1 ,x 2 ,t) = 



1 



T-t 



tfj(x 1 ,t r )ij(x 2 ,t' +t) dt', (3) 



where T is the duration of the observation and if> is the ob- 
served signal. Pre-multiplication by the temporal window 
function (equal to unity in the interval [— T/2, T/2] and zero 
outside) is included in the definition of tp. The positive time 
lags (t > 0) give information about waves propagating from 
X\ to X2, and the negative time lags (t < 0) give informa- 
tion about waves propagating in the opposite direction. The 
cross-covariance function provides a means to add random 
waves constructively; it represents, in some sense, a solar 
seismogram. 

The goal is to find a three-dimensional model of the 
Sun which is consistent with a large set of observed cross- 
covariance functions. Full waveform modelling has not 
been attempted yet. In practice, only a few parameters are 
used to describe a cross-covariance function. The most im- 
portant one is perhaps the travel time. 
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Two methods are used to measure travel times from 
C(t). The first method [(Duvall et al. 1997)1 consists of fit- 
ting a Gaussian wavelet to each branch of the cross- 
covariance. For example, the positive-time part of the cross- 
covariance is fitted with a function of the form 



samples of the data (Jen sen, Duvall & Jacobsen 2003)| We 
note that the noise covariance matrix depends on the fourth 
moments of the wavefield, 



Aexp[-7 2 (t - i g ) 2 ] cos[uj (t - 7+)] 



(4) 



where all parameters are free. Some spatial averaging is of- 
ten necessary before the fit is performed to guarantee con- 
vergence. The time t + is called the phase travel time. Simi- 
larly a phase travel time r_ is measured by fitting a wavelet 
to the negative-time part of the cross-covariance. The travel- 
time differences ram — t + — t_ are sensitive to internal 
flows, while the mean travel times T moan = (r+ + t_ ) /2 are 
sensitive to wave-speed perturbations (temperature, density, 
magnetic field). This distinction, however, may be fuzzy due 
to finite-wavelength effects. Note that the phase travel time 
is defined modulo 2it/ujq. 

An alternative definition of travel time was introduced 
by Gizon & Birch (2002, 2004) by analogy with a defini- 
tion used in geophysics (Zhao & Jorda n~T998)[ This method 
is designed to measure travel times from cross-covariances 
measured with short T and with as little spatial averaging as 
possible. By definition, the travel times r + and t_ are given 
by 

r±= Jw ± (A,t)[C(x 1 ,x 2 ,t)-C (A,t)]&t, (5) 

where the function Cq is a smooth reference cross- 
correlation derived from a horizontally homogeneous solar 
model, A = Xi — x\, and the weight functions W± are 
proportional to <9 t Cq. One advantage of this definition is a 
straightforward, linear relationship between travel times and 
cross-covariance functions. 

The second definition of travel time (eq. O) is use- 
ful as it is more robust to noise than the wavelet fits 
(Rot h, Gizon & Beck 2007)] It is also much easier to 
interpret in terms of perturbations to a reference solar 
model (sec 3]). On the other hand, a wavelet fit returns 
additional information: the amplitude, A, the group 
travel time, t s , the frequency, ojq, and a measure of 
the width of the wavelet, I/7. Amplitudes, in particu- 
lar, are known to vary very significantly in the vicinity 
(Jensen, Pijpers & Thompson 2006) 



of active regions 

since sunspots 'absorb' incoming acoustic energy 



(Braun, Duvall & Labonte T987)] |Woodard (2002)] and 
Rajaguru et al. (2006a) provide two different reasons 



why amplitude variations and travel-time shifts are not 
independent parameters as a result of inhomogeneous wave 
damping and excitation. 

The main source of random noise in local helioseismol- 
ogy is, by far, the stochastic nature of solar oscillations. A 
good understanding of the properties of noise is needed for 
a correct interpretation of the measurements and, in partic- 
ular, for solving the inverse problem (see sect. 0. Noise is 
specified by the covariance matrix of the travel times, which 
may be estimated directly by spatial averaging over many 



Cov[<£*(fci,a;i)(/>(fc2, w 2 ),</>*(fe3,w 3 )<Xfe4, w 4 )]. 



(6) 



It may be written down explicitly (Gizon & Birch 2004) 
when the cf>(k,oj) are assumed to be normally distributed 
and uncorrelated in Fourier space, i.e. under the assumption 
that the random function <p(x, t) is spatially homogeneous. 
This last assumption, however, is not valid near active re- 
gions: noise in travel-time measurements also depends on 
the local solar properties that we are trying to infer from the 
travel times (this has been happily ignored). 

To summarize this section, the challenge is to extract 
unbiased information from the wave field, which is relevant 
and can be interpreted. This can hardly be achieved with- 
out some understanding of the physics of wave propagation 
through complex media, which we now discuss. 



4 Forward modelling 

Forward modeling is the process of computing the wave 
field given a prescribed solar model. The solar interior is 
a complex medium in at least two different ways. The no- 
tion of complexity depends somewhat on which waves we 
are talking about. First, the solar convection zone has in- 
homogeneities that are smaller than the wavelengths of so- 
lar oscillations and/or that vary faster than the wave peri- 
ods. For such a medium, ray acoustics is not valid. The sec- 
ond kind of complexity is that some inhomogeneities, like 
sunspots, are not small-amplitude perturbations to the aver- 
age medium, so that single-scattering approximations can- 
not be used. In the case of turbulent convection, both levels 
of complexity are combined. 

Needless to say that a lot of work remains to be done in 
order to model wave propagation through complex/random 
media in the context of local helioseismology. For the sake 
of simplicity, it is often assumed that there exists a sepa- 
ration of scales between turbulent convection, waves, and 
slowly-varying structures. The effects of turbulent flows are 
twofold: excitation and damping of the waves. Both effects 
are usually parametrized using phenomenological models. 
Nearly all studies have assumed that the solar waves (~ 

5 min periods) propagate in a background that is steady 
over the duration, T, of the observations (typically several 
hours). 

So far, most efforts have focussed on wave propagation 
through a steady, weakly inhomogeneous background spec- 
ified by small-amplitude perturbations to a reference solar 
model invariant under horizontal translation. This class of 
problems has been studied extensively in geoseismology. 
The basic idea is to use first-order perturbation theory - the 
first Born or Rytov approximations - to write the scattered 
wave field in terms of the zero-order wave field. 
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The equations of solar oscillations (see e.g. Cameron, 
Gizon & Daiffallah 2007) may be combined in the form 
|(Lynden-Bell & Ostriker 1967) 

m +m = s, (7) 

where L = pdf + ■ ■ ■ is a linear differential operator 
governing adiabatic oscillations, T is a phenomenological 
damping operator, and S is a stochastic source function 
(granulation). Using subscripts to label the unperturbed 
quantities defined in the reference solar model, we have 
(Lq + To)[£ ] = So. The background solar model is as- 
sumed to have no magnetic field and no flow. The random 
source function, So, is assumed to be horizontally homo- 
geneous and stationary in time. For a background model 
where density and pressure only depend on height, £ can 
relatively easily be written in terms of a Green's tensor us- 
ing a normal mode summation (Birch, Kosovichev & Duvall 
2004). In the first-order Born approximation, the first order 
perturbation to the wave displacement, <5£, is such that 

(L + T )[6t] = -(6L + 6r)[S ]+6S. (8) 

Inhomogeneities, for example sound speed variations 6c(r) 
or a buried magnetic field B(r), are encapsulated in the 
wave operator SL. Note that we have included the pertur- 
bations to the damping operator and to the source function. 
These two terms should not be neglected a priori. For exam- 
ple, it is not excluded that a 5S could arise from a change 
in the properties of granulation due to a local change in 
sound-speed (or temperature). These terms, though, are un- 
likely to be dominant, except perhaps in magnetic regions 
(see previous section). The unperturbed Green's tensor can 
be used to solve equation ([8]) for 6£, just like it was used to 
solve the zero-order problem. The Born approximation is an 
equivalent-source description of wave interaction (the right 
hand side of eq. [HI). It has been used to compute the effect 
of sound speed perturbations (Birch, Kosovichev & Duvall 
2004) with 

- 5L[£ ] = V(p 5c 2 V ■ £ ). (9) 

In this volume Birc h & Gizon (2007)| consider, for the first 
time, the effect of a flow, U(r), on wave travel times using 

- 6L[£ ] ~ 2 Po U ■ Vd t £ . (10) 

In the case of magnetic perturbations, the appropriate oper- 
ator to consider (Lorentz force) is 
1 



8L[t ] = -^(J'*B- 



J x B') 



(11) 



where currents are given by J = V x B and J' — V x B' 
and the magnetic field oscillations are B' = V x (£ x B). 
Note that to first order the perturbation is quadratic in the 
magnetic field. 

Once the wave displacement, £ = £ + S£, has been 
computed, it is possible to obtain the first-order perturbation 
to the cross-covariance, SC = C — Co, using equations (Q~|i- 
(0). The next step is to relate the observations, in particular 
the travel times, to SC. This operation is trivial with the one- 
parameter fit since St = J W 8C dt. It is much less obvious 
with the (non-linear) wavelet fit. 



The result of single-scattering theories is that there ex- 
ists a linear relationship between travel-time perturbations 
and small-amplitude changes in internal solar properties. It 
is helpful to cast the observational constraints into the form 



8n = 



K?(r)5q a (r) dr, 



(12) 



where i indexes the M different measured travel times (so i 
runs from 1 to M). The index a refers to the various types 
of physical quantities, q a , like sound speed, density, first 
adiabatic exponent, flows, magnetic field, etc. The functions 
Kf(r) give the sensitivity of small changes in travel time, 
Sri, due to a change Sq a (r) throughout the solar volume. 

The accuracy of the Born approximation 
for magnetic perturbations has been tested by 
Gizon, Hanasoge & Birch (2006) using the exact so- 



lution for waves impacting a magnetic cylinder in an 
otherwise uniform medium. For a one kilogauss magnetic 
field, the Born approximation would appear to be valid 
except close to the solar surface (the first few 100 km). 
The assumption of small perturbations breaks down in 
active region sub-photospheres. At this point a treatment 
of the strong perturbation regime is needed, a major topic 
of current research. A very good account of the types of 
problems encountered in the magnetoseismology of active 
regions is provided in this volume by Cally (2007) using an 



extension of ray theory to magnetoacoustic waves. We refer 
the reader to this paper and references therein. 

5 Inversions 

The linear inverse problem consists of making inferences 
about the unknown functions Sq a (r) from the measure- 
ments <5tj. We assume (see previous section) that the kernel 
functions Kf(r) are known sufficiently accurately. 

Equation (TT21) can be written in a vector notation, with a 
dot product defined in an obvious way, as 



K • Sqdr 



or equivalently 



Ki • (5qdr. 



(13) 



(14) 



In practice, all the data analysis is done on computers 
where quantities are represented in discretized form. Even 
if that were not so, various methods of solving the inverse 
problem involve making a discretized representation over a 
finite basis of functions. Hence we write 



d = Ai 



or 



(15) 



(16) 



where di represent the data, Xj is a discrete representation of 
the unknown functions, and multiplication by Aij replaces 
the integral with weighting given by the kernel functions. 
Let us say that j runs from 1 to N. 
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Taking a linear combination of the data constraints, 
weighted by some as yet arbitrary coefficients Cj, one ob- 
tains 



Cidi = 




C,;K ?; 



<5q dr. 



(17) 



The combination K, = CjKj is called the averaging ker- 
nel. Typically one attempts to localize 7C near some target 
location ro; in the case that there are several components a, 
localization means localizing the component corresponding 
to one chosen a whilst simultaneously making the compo- 
nents of 7C corresponding to all other as small. There are 
various ways to choose the coefficients Ci, different ways 
leading to different 'solutions' of the inverse problem. The 
set of such methods forms a family of optimally localized 
averaging (OLA) methods. An OLA procedure for the in- 
version of f-mode time-distance helioseismology is given 
by Jackiewicz et al. (2007). 

A very different approach, at least philosophically, is 
'data fitting', by finding the solution Xi to minimize the data 
mismatch 



Id- Axil 



(18) 



In practice, strictly minimizing the data mismatch typically 
leads to an unacceptable solution x for noisy data, because 
matrix A is usually ill-conditioned. The problem must be 
regularized in some way. One way that is common in global 
helioseismology but less so in local helioseismology is so- 
called Tikhonov regularization, in which one minimizes not 
simply the data mismatch but the sum of this plus a penalty 
term that gets large as x exhibits some undesirable prop- 
erty such as having a large magnitude or large fluctuations. 
Another way is to solve the minimization problem ( TT8l in a 
succession of iterative steps that would eventually converge 
to the exact solution, but instead truncate the process after 
a small number of iterations before the exact solution has 
been reached and, more importantly, before the solution in 
attempting to fit closely the noisy data has developed un- 
acceptable features. Such a method is the so-called LSQR 



(Paige & Saunders 1982) 



Complementary to the inversion methods above is so- 
called forward-modeling, in which one searches through 
a space of possible models to find a best fit to the ob- 
servational data. One such search method is genetic mod- 
eling, which has been used in global helioseismology by 
Char bonneau et al. (1998) and recently in local helioseis- 
mology of sunspots by |Crouch et al. (1 990) 

One consideration in picking an inversion method is 
how long the solution will take to compute, or indeed 
whether the problem will even fit within the capacity of the 
available computers. Typically, solving an OLA problem in- 
volves inverting one or more matrices of size M x M. The 
data-fitting techniques involve inverting one or more matri- 
ces of size M x N. If M and N are very large but one 
is much larger than the other, the difference between these 
two matrix inversion problems may be a significant practi- 
cal consideration. 



Time-distance helioseismology typically involves 
taking measurements of travel times over a variety of 
distances, at many locations on the visible disk of the Sun. 
If the kernels are invariant under horizontal translation 
of the measurement point, then this can be exploited 
to speed up the inversion. A way of doing this is mul- 
tichannel deconvolution (MCD) which was introduced 
into helioseismology by B. H. Jacobsen and J. M. Jensen 
(Jensen, Jacobsen & Christensen-Dalsgaard 1998) After 
Fourier-transforming the problem in the horizontal, this 
reduces the problem to a series of one-dimensional inverse 
problems in the depth direction. A refinement to control 
regularization in the horizontal direction has been explored 



by|Zharkov, Nicholas & Thompson (2006) 



Unfortunately, whilst it is approximately the case that 
the problem is invariant under horizontal translation near 
the centre of the disk, foreshortening and line-of-sight 
effects mean that this assumption breaks down as one 
moves away from disk centre. Kernels for time-distance 
helioseismology using f-mode data have been calculated 
by [Jackiewicz , Gizon & Birch (20 06) demonstrating ex- 
plicitly that these vary with horizontal location on the disk. 
The kernels can be expressed as a weighted sum of a num- 
ber of different contributions: each of the individual contri- 
butions is the same at all locations but the weight given to 
each depends on location. Whether the computational effi- 
ciencies of MCD can somehow be retained when the kernels 
are no longer translationally invariant is an open question. 

A similar concern arises if the initial background model 
is horizontally inhomogeneous. This would be the case, for 
example, if one were to proceed with the inversion itera- 
tively, updating the background model to incorporate the 
current inversion results after each linear inversion step. 

In common with global helioseismic studies, an is- 
sue in local helioseismology is how to test the reliabil- 
ity of the inferences drawn from inversions. The inver- 
sion solution depends on the method chosen and the val- 
ues selected for various trade-off parameters. These af- 
fect the resolution attained and how errors propagate from 
the data to the solution. For a strictly linear problem and 
linear inversion method, the resolution can be understood 
completely in terms of averaging kernels: the solution at 
each point can be rigorously interpreted as a convolu- 
tion of the true underlying solution with the averaging 



kernel (Christensen-Dalsgaard, Schou & Thompson 1990) 
Resolution measures can be defined based on properties of 
the averaging kernels. Also, for a linear inversion method, if 
the statistics of the data errors are completely known, then 
the statistics of the errors on the solution can be determined 
precisely, including the correlation between errors at differ- 
ent points in the solution. Complicating factors that are less 
easy to account for are: systematic errors in the data; inade- 
quacies in the description of the data errors; nonlinearity of 
the underlying problem; inaccuracies of the forward model 
(as expressed in the kernels), which point could include the 
issue of nonlinearity. The effects of some of these might be 
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tested by making inversions of artificial data generated from 
more or less sophisticated numerical simulations of waves 
in models of the whole or part of the Sun. 

Another issue is that many different physical effects en- 
ter (as indexed by a in eq. |[T2"1 ). To what extent can these 
be separated by the inversion, either in principle, or in prac- 
tice given a finite number of noisy data? As discussed by 



Zweibel & Gough (1995) in the global helioseismic con- 
text, it is not possible for example to distinguish unambigu- 
ously between thermal perturbations and magnetic fields 
with a finite number of frequency measurements. This is an- 
other aspect of the ill-posedness of the inversion problem. 
Other a priori knowledge, or assumptions or prejudices, 
might be introduced to break the ambiguity. It may also 
be possible to supplement the seismic data with other ob- 
servational constraints, such as photospheric magnetic field 
measurements for instance. 

6 Discussion 

We have tried to provide a short description of some of 
the current topics of research in time-distance helioseismol- 
ogy. Although much work remains to be done, the general 
procedure to be followed to infer small-amplitude, steady 
perturbations in the solar atmosphere is, in principle, rel- 
atively well understood. The treatment of strong, complex 
perturbations is, however, still in its infancy. As discussed 
by Werne, Birch & Julien (2004), the study of wave prop- 
agation through magnetic active regions would benefit a 
lot from numerical simulations. Indeed large-box realis- 
tic simulations of fully-compressible non-linear magneto- 
convection are becoming feasible (e.g. Georgobiani et al. 
2006; Steiner 2007) offering test beds to validate the var- 
ious methods of local helioseismology (Zhao et al. 2006). 
Other tests will be performed using codes that propagate lin- 
ear waves through inhomogeneous background solar mod- 
els (e.g. Cameron, Gizon & Daiffallah 2007; Hanasoge & 
Duvall 2007), which are much less computer intensive. 

Much progress in local helioseismology is to be ex- 
pected in the next few years thanks to numerical mod- 
elling of wave propagation together with insight gained 
from asymptotic methods (e.g. Cally 2007; Gordovskyy & 
Jain 2007; Gough 2007) and from other fields of physics 
(e.g. seismology, ocean acoustics). Refined methods of data 
analysis, such as e.g. the three-point correlation technique 



of Pijpers (2007) may also prove very useful in extracting 
all the information from the MDI/SOHO and the upcoming 
HMI/SDO observations. 
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